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Transition path sampling is a rare-event method that estimates state-to-state time- 
correlation functions in many-body systems from samples of short trajectories. In 
this framework, it is proposed to bias the importance function using the lowest Jaco- 
bian eigenvalue moduli along the dynamical trajectory. A lowest eigenvalue modulus 
is related to the lowest eigenvalue of the Hessian matrix and is evaluated here using 
the Lanczos algorithm as in activation-relaxation techniques. This results in favor- 
ing the sampling of activated trajectories and enhancing the occurrence of the rare 
reactive trajectories of interest, those corresponding to transitions between locally 
stable states. Estimating the time-correlation functions involves unbiasing the sam- 
ple of simulated trajectories which is done using the multi-state Bennett acceptance 
ratio (MBAR) method. To assess the performance of our procedure, we compute 
the time-correlation function associated with the migration of a vacancy in a-iron. 
The derivative of the estimated time-correlation function yields a migration rate in 
agreement with the one given by transition state theory. Besides, we show that the 
information relative to rejected trajectories can be recycled within MBAR, resulting 
in a substantial speed-up. Unlike original transition path-sampling, our approach 
does not require computing the reversible work to confine the trajectory endpoints 
to a reactive state. 
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I. INTRODUCTION 

Rate constants of thermally activated processes are crucial parameters governing kinetic 
and transport properties of condensed matter. Estimating these rates directly using molec- 
ular dynamics or kinetic Monte Carlo typically remains a challenge as the probability of 
observing an event is very low at the simulation time-scale. In practice, a sufficient number 
of reactive events has to be collected so as to achieve high accuracy in the estimates of 
the rates. This task consumes considerable amounts of computational time or may even 
be prohibitive. This limitation arises in the simulations of many physical, chemical or bio- 
logical systems in which the structural evolution of clusters, molecules or proteins involves 
transitions between locally stable states occurring rarely during a computer simulation. To 
reduce the total simulation cost, several biased sampling methods have been developed in 
order to enhance the probability of observing the rare events.-"- Among them, transition 
path sampling^ (TPS) aims at estimating a time-correlation function between two indicator 
functions, Ha and hs-, characterizing the reaction. The indicator function /ia(_b)(x) takes 
value 1 if state x G A{B) and value elsewhere, A and B being two basins traditionally 
referred to as reactant and product. In practice, the state-to-state time-correlation that is 
estimated in TPS is based on a path ensemble average rather than a time average. De- 
noting hy z = {x(s)}o<s<r the trajectories of the path ensemble and of duration T, TPS 
considers that the trajectory initial states x(0) are distributed according to Z^-aP^' where 
denotes the equilibrium Boltzmann distribution normalized to one within basin A. The 
time-correlation function in TPS is the overall conditional probability that a trajectory from 
basin A is in state B aX t < T 



where Pcond(^^-2|x(0)), the probability measure of trajectory z, is conditioned on its initial 
state x(0). The probability measure of the initial states is [x(0)] [(ix(0)]. This defi- 
nition of the time-correlation function^»^ naturally applies to any stochastic model coupling 
the dynamics to a thermostat, such as the Langevin dynamics. It may also apply to a model 
in which the system is initially prepared in equilibrium and allowed to evolve in isolation 
(without exchanging heat with the environment) by a deterministic reversible dynamics. 
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The latter kind of dynamical model is often studied as it well approximates situations where 
the times t of relevance to rate processes are small compared with the time required for a 
significant heat exchange of the simulated system with the environment. 

The conditional probability C{t) plays an important role in the molecular simulations of 
rare events because it corresponds to the fraction of trajectories that are reactive. When 
fast molecular relaxations within basin A and activated processes associated with transitions 
exiting basin A occur on well-separated time scales, the derivative of the time-correlation 
function, dC (t) / dt, displays a plateau corresponding to the phenomenological rate constant 
kA^B- The phenomenological rate may be derived at macroscopic scale considering fluxes 
over populations or at atomic scale resorting to linear response theory. TPS estimates 
the time-correlation function C{t) over a time interval ]0,T] from samples of constrained 
trajectories. In practice (see Refs. LlJ and la), TPS factorizes C{t) into the product of (i) 
a reduced time-correlation function Rt>{t) = C{t)/C{t'), estimated by taking an average in 
the ensemble of transition paths and (ii) a constant factor Kf = C{t'), obtained from the 
reversible work required to progressively confine the trajectory endpoints beyond milestones 
marked out along the reaction up to the reactive state. 

We show herein that, for a thermally activated process, it is possible to compute the 
time- correlation function C{t) without confining the trajectory endpoint to a reactive state. 
Eigenvalues of the Jacobian matrix associated with the dynamics will be used to favor the 
sampling of active trajectories, those visiting negatively curved portions of the potential en- 
ergy surface (i.e. saddle regions of instability between energy basins). This strategy involves 
post-processing the samples of active trajectories so as to simultaneously (i) correct for the 
bias and (ii) estimate the fractions of reactive segments yielding the entire correlation func- 
tion C{t), two requirements reflecting the two-step procedure of TPS. We show in the paper 
that the multi-state Bennett acceptance ratio method^ achieves this goal. The proposed 
method is referred to as SUNDAE as it consists in sampling and unbiasing dynamically 
activated events. 

This paper is organized as follows. The key ingredients of the methodological recipe 
are given in the five following sections. Section |TT] derives the conditional probabilities for 
generating trial trajectories. Section [TTll describes the Verlet map. Section [IV] decomposes 



the Jacobian matrix associated with the Verlet map. Section |V] describes the path-samphng 
scheme and section |VT] derives the various unbiasing estimators tested in this study. We show 
in particular that MBAR can be combined with waste-recychng,— 1^2, a technique consisting 
in including information from unselected proposals so as to reduce the statistical variance of 
the estimates.— SUNDAE is finally applied to the calculation of a rate constant associated 
with the migration of a vacancy in a-iron (Sec. IVIip . 

II. TRAJECTORY SPACE AND CONDITIONAL PROBABILITIES 

Let {qi}i<i<3N and {pi}i<i<3N denote the canonical coordinates of the system. In the 
following we write the evolution equations in matricial form and denote vectors and ma- 
trices by lowercase and uppercase bold letters, respectively. We thus introduce the column 
vectors q and p for the positions and momenta. State x is also written as a column vector 
(gi, g3Ar,pi, ...,p3Ar)^, the upper T denoting vector or matrix transposes. Vector m de- 
notes a deterministic and invertible map defined on E, the phase space: x G -E — t- m(x) G E. 
The inverse map is denoted by m^^. It is also a functional column vector. The Jacobian 
matrix of m is denoted by (Vm)-^ : its components are [(Vm)^]^_^. = Vj-rrij where Vj rep- 
resent the partial derivative From the inverse function theorem, the Jacobian matrix of 



Consequently, the Jacobian determinant of m ^ at m(x) is det[Vm] ^ at x. Let Sa{dy) 
denote the Dirac measure at a. The direct transition probability from x to y is 



Similarly for the inverted map, the probability to transition from y to x (reverse transition) 
is 



where the substitution of m~^(y) by x in the determinant does not affect the distribution. 

To explain the physical meaning of the Jacobian determinants, let consider a small hy- 
perparallelepiped represented by one corner xq and a matrix 6\o wherein the j-th column 




6^ (m-i(y)) = |det [Vm(x)]| (y) . 



(2) 



Sy [m(x)] = |det [Vm(x)]| 5m-i(y) W 



(3) 
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vector, denoted by 5xq is the j-th edge arising from xq. By successive applications of the 
map, the j-th vector evolves according to 5x:^_,_]^ = m(xn + 6x.D — m(x„) where 5x^ and 
S^n^i are the j-th edge vector at x„ and x„+i = m(x„). Expanding the map m at x„ 
yields^ 

m(x„ + 5<) = m(x„) + (Vm)'^(x„)5x^„ + 0{\\6^if). (4) 

Thus the j-th column vector evolves according to ^x^^^^ = Vm^(x„)5x{ at first order. The 
matrix representing the small hyperparallelepiped similarly evolves according to S\n+i = 
Vm"^(x„)5V„. Now, if matrix components at step are chosen such that det(5Vo) > 0, 
then this quantity characterizes the initial volume of the small hyperparallelepiped. This 
volume then evolves according to 

det(5V,) = det(5Vo) Ui=\ det [Vm(x„)] . (5) 

Thus, the Jacobian determinant product represents the factor by which phase space is com- 
pressed or expanded on the £ successive applications of the map m starting at xq. 

A trajectory is defined as an indexed sequence of L + 1 states in phase space, denoted 
hj z = {xq, xj^}, independently from the map m. From (|2]) and (|3]), the conditional 
probability Pcond(-2|x^) to generate z by successively applying the map starting from state 
X£ (0 < £ < L) is the product of the forward and backward transition probabilities from the 
corresponding index i 

Pcond(2;|x^) = Pf^d{z\i, yie)Phwd{z\i, x^.); (6) 
Fi„diz\i,Xi) = ni=l |det [Vm(x„)]| 5rn(x,0 (Xn+l) , 
Pbwd(^|^,x^) = nl=l) |det [Vm(x„)]r^ '^m-Hx^+i) M ■ 

For Dirac measures, we have equivalence between 6a{dh)da and 6h{daL)dh. Hence, the 
conditional probability (|6]) can be cast in the general form 

Pcond(2;|x^) = Pcond(^|xo) 111=0 [G(x„)]"^ , (7) 

where G = (Vm)^Vm denotes the Gram matrix. This identity is a key ingredient of the 
path-samphng schemes described in Sec |Vl We now focus on the map that will be used. 
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III. VERLET MAP 



We consider djTiamics governed by an Hamiltonian 'H(x) = /C(p) + V(q) where V(q) 



and /C(p) = \Y^i=iP'i/^i denote the potential and kinetic energies, respectively. The 
Hamiltonian gradient V'H(x) is written as a column vector (ViH, VeArH)^. Hamilton's 
equation expressed in matricial form reads^ 

X = JV?{(x) (8) 

where J is the canonical skewed- symmetric matrix 




and I3JV is the 3A^ x 3A^ identity matrix. Hamilton's evolution equation is discretized through 
a map hIt- by approximating x(rir) with x„ = mT-(x„_i) where r is the time-step. Trajec- 
tories are z = {xo,...,X2,} = {x(0), x(T)} where the total length T = Lr. Let define 
the vectorial functions q and p by q(x) = q and p(x) = p and let o denote the function 
composition. We write /C = /C o p and V = V o q, making it possible to apply the general- 
ized gradient V on the kinetic and potential energies directly. The position- Verlet map is 
obtained by discretizing Hamilton's equation as follows 

mr= (i+ i JV^) 0(^1 + r JVV) o (^i + f JV^) . (9) 

where i is the functional identity : i(x) = x. Defining qn+1/2 = q o (i + | JV/C) (x„) 
enables one to decompose x„+i = m^(x„) according to the intermediate updates of the 
splitting (ED 

qn+1/2 = Qn +iVp/C (p„) , (10a) 

Pn+l = Pn -rVqV(q„+l/2), (10b) 

qn+1 = qn+1/2 +iVp/C (p„+i) , (lOc) 

where qj,n+i simplifies into g^-n + ^{pj,n + Pj,n+i)/mj. 
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IV. EIGENVALUE DECOMPOSITION 



Verlet maps have three interesting properties'^ : (i) time symmetry (i.e. = m~'); 

(ii) reversibihty [i.e. = r o 111,- o r where r is the momentum reversal apphcation: 
i^yP^)'^ — i- (q-'", — P"^)"^] and (iii) symplecticity [i.e (VmT-)"^JVmT- = J]. To prove property 

(iii) , we formally expand the Jacobian matrix for the position- Verlet map 



where K and V are the Hessian matrices of /C, V, and I = Vi^. The components are 
(K.)ij = VjVjTC, (V)jj = VjVjV and = 6ij. In (ITT]) the Jacobian matrix (Vm^)-^ 
is factored into a product of three symplectic matrices, hence it is also symplectic. Note 
that the Hessian of the potentiel energy is to be evaluated at qn+1/2 for (Vm^)^(x„) or 
(Vm_^)-^(x„+i) as the potential gradient in f lTOj) . 

An important implication of the symplecticity property is that eigenvalues of the Jacobian 
matrix occur in reciprocal pairs: if fi is an eigenvalue, then so is fi~^ with the same algebraic 
multiplicity. To facilitate the eigenvalue decomposition in reciprocal pairs, we perform a 
change of basis in order to partition the Jacobian matrix into four commuting and symmetric 
block matrices. The associated transformation matrix is defined by 



The transformation amounts to switching to mass-weighted coordinates. The two block 
matrices are diagonal and have dimension 3N x 3N with entries (K")^^- = m~°'Sij. Note 
that K' is the Hessian matrix of the kinetic energy /C. After expressing the various matrices 
composing the Jacobian matrix (Eq. flTT]) ) in a basis with mass- weighted coordinates, the 
latter one becomes r^'(VmT-)-^r and writes explicitly 



(Vm,)^ = (I + iJK)(I + rJV)(I + f JK), 



(11) 





(12) 



2 1 1 

with C = IsAT — VK2. The mass-weighted Hessian matrix of the potential energy V is 
symmetric since we have 
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Hence, Kz VK2 admits real eigenvalues, denoted by Uj and sorted here in ascending order. 
An unstable normal mode is characterized by a negative eigenvalue u'^ in which case the 
eigenfrequency uj is pure imaginary. A stable normal mode is characterized by a positive 
eigenvalue Uj in which case the eigenfrequency Uj is real. The eigenvalues of C are Cj = 
1 - u]ry2. 

The commutative subring formed by the block matrices C, f (C + Isa?) and ^(C — I^n) 
within the partitioned Jacobian matrix (fT2|) allows us to factorize its characteristic polyno- 
mial (secular equation) as follows (/x G C) 

det [(/ilsiv - C)' + 13^ - C^] = 0. (13) 



Diagonalizing the matrices of f[T3|) in the real eigenbasis of C enables one to decompose the 
secular equation into the product of the following second order equations 

(/i - Cj)2 + l-c'j=0. (14) 

If fXj is a root of ( 11^ . then Vieta's formula implies that is the second root and we have 



For convenience, we introduce an effective frequency uj to express the matching pairs of 
roots as 

= exp {±iujjT) . 

For Cj > 1 (i.e. Uj < 0), uj and uj are both pure imaginary numbers. We have 



cosh(iujjT) = Cj and smh(icojT) = — ^ / |1 — < 0. Consequently, icoj < and smh.{iujjT/2) 



— y ^[cosh.{iojjT) — 1] = ~\J\{cj — 1). Substituting 1 — a;|r^/2 for Cj in the last square root 
and setting uoj to the value i\ujj\ yields smh{iujjT /2) = iujjT/2. The corresponding Jacobian 
eigenvalues jjij and are real. The projections of 5x along the eigendirection of jjij and 
are contracting and expanding, respectively. The occurrence of an unstable normal mode 
characterizes configurations in saddle regions. For Cj = 1, we have uj = uj = 0. In any 
energy surface exhibiting periodic boundary conditions (as will be the case in Sec. IVIip there 



is a zero frequency for each normal mode associated with a translational symmetry. For 

8 



-1 < Cj < 1, (i.e. < < 4/r^), Uj and uj are both real positive. We have cos(a;jr) = Cj 



and sin(wjr) = yl — > 0. Since uj is chosen positive, we have sin(a;jr/2) = UjT/2. 
In this case, is a conjugate pair in the unitary circle of the complex plane (with 

unit modulus). This favors oscillations of the projection of 5x into the plane containing the 
two corresponding eigendirections. States in metastable basins have all their normal modes 
stable {u'j positive). For Cj < —1 (i.e. > 4r~^), coj is pure imaginary while the corre- 
sponding normal mode is stable {uj is real positive). The effective mode of the discretized 
dynamics is unstable but becomes stable again if the time-step r is decreased below 2/uj. 
The situation cj < —1 thus involves numerical instability and will not be covered in the 
following as admissible time-steps in molecular dynamics should be a fraction of the fastest 
vibration period (see the linear stability analysis in Ref. 
Hence, assuming w| G M and w| < 4r^^, we may choose 
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UJj = — 

T 



arcsin oRe (^^j'^^ + i x arsinh o Im (^^j'^^ ; (15) 



where Re and Im denote the real and imaginary parts of a complex number and o still denotes 
the functional composition. Equation (fT6!l generalizes in the complex plane as follows 



UJj = 2(zr)-^ In (^^I-u^t^A + iUjT/2^ , (16) 

provided the principal branch of the multi-form complex function is delimited by the cutting 
segments ] — oo, — 2r~^[ and ]2r~^, +oo[. The linear approximation of the effective frequency 
Uj is Uj. To detail the higher terms, let first derivate f lT6|) with respect to Uj and then expand 
the result in Taylor series (wj G M and Uj < 4r^^): 



l-uy/A = J20{\u,rY'. (17) 



da;,- / . 

k=0 



The effective frequency Uj is obtained by integrating the function u — )■ J2T=o i^k) i^'^'^Y'^ 



from to Uj 



fc=0 

where we used the fact that Uj and Uj both take value zero at Cj = 1. The series is real 



positive whatever cj^ G M and cj^ < 4r ^. The non-linear contributions absorbed into u 



contain even powers of r exclusively, a property reflecting the oddity of the circular and 
hyperbolic inverse sines in f|T5|) . The first terms of the series are 



Uj = UJj 



' + 24^^- + 640^^- + 7168^^- + 



In practice, an appropriate truncation of (fT8|) is used to avoid numerical round-off errors 

when evaluating inverse sine functions close to 0. 

Let us now characterize the extremal moduli of the eigenspectrum. The imaginary parts 

Im(ci;j) form a descending series, as for the imaginary parts Im (wj) since duj/dujj is real 

positive in (ITTj) . Conversely, the series formed by the moduli = exp oIm(— cD^r) is 

ascending for the exponential function being increasing on M with lm(— cDj). Therefore, the 

smallest and largest moduli are and two quantities always associated with the 

lowest eigenvalue cu^. Physically, the matching pairs of Jacobian eigenvalues are 

local Lyapunov numbers. The local Lyapunov exponent, defined by 

~ 2 r 
— r In = Im(a;i) = — arsinh o Im(— tui) 

is associated to the maximum expansion rate of a perturbation 5x„ at current state x„. 
Lyapunov exponents rather correspond to the expansion rate in the infinite time limit r 

Lyapunov exponents measure the sensitivity of dynamical systems to small changes in 
initial conditions and, for this reason, are widely employed for characterizing the dynamical 
properties of many nonlinear and Hamiltonian systems.— For instance, by monitoring the 
Jacobian eigenvalues of the velocity- Verlet map and averaging the logarithm of the eigenvalue 
moduli exceeding one,—"— the chaotic or regular nature of dynamical trajectories could be 
identified and related to the occurrence or absence of structural transitions in small Lennard- 
Jones clusters. In two recent investigations, the Jacobian eigenvalue spectrum itself has 
been successfully used to enhance the occurrence of rare chaotic and regular trajectories in 
computer simulations of dynamical systems.— In the Lyapunov-weighted dynamics (LWD) 
developed by Tailleur and Kurchan^ (see also Refs. |24| and picciani2011simulating), a set of 
small vectors evolves according to the Jacobian equation (jl]) coupled to a birth-death process 
consisting in multiplicating or deleting vectors depending on their norms. In the Lyapunov- 
weighted path-sampling (LWPS) method developed by Geiger and Dellago,'2S the birth-death 
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process is replaced by a Metropolis algorithm sampling short dynamical trajectories. The 
Metropolis acceptance rate of LWPS includes a biasing factor proportional to the relative 
Lyapunov indicator^^ of the involved trajectories. The ability of LWPS to identify rare 
trajectories suggests that this approach is well suited to the estimation of time-correlation 
functions associated with thermally activated processes. 



V. SAMPLING DYNAMICALLY ACTIVATED EVENTS 

We implement a variant form of LWPS. Instead of resorting to the relative Lyapunov 
indicator we use Y[n=o lA*i,n+i/2|, the product of the smallest eigenvalue moduli along tra- 
jectory z, to bias the importance function. The lowest eigenvalue /ii,n+i/2 refers to Jacobian 
matrices (VmT-)-^(x„) or (Vm_T-)'^(x„+i). The logarithm of the biasing quantity reads 



C{z) = -2 ^ arsinh o Im (^a;i(q„+i/2)0 < 0. 



n=0 



C{z) is hereafter called "activation indicator". By definition, a trajectory is said to be 
inactive when its indicator is zero and to be active when it is strictly negative, a negatively 
curved portion of the energy surface being visited in the latter case. At variance, the relative 
Lyapunov indicator does not always discriminate active trajectories from inactive ones. 

The lowest eigenvalue uf (qn+1/2) is calculated using the Lanczos algorithnt^i as in 
activation-relaxation techniques.—"— At each implementation, the algorithm finds the lowest 
(or highest) eigenvalue and eigenvector of any matrix at a reduced computational cost, by 
constructing a Krylov basis of size d and diagonalizing the associated d x d matrix where 



33 



d is small (see for example appendix A of ReL^- for details). As pointed out in Ref. 
using a Krylov basis of size d = 15 guarantees the convergence of eigenvalue cof (qn+i/2)- 
Herein, we have decreased the basis size to d = 8 to increase numerical efficiency. At each 
implementation, we verify that the Lanczos solution is stable; if it is not the case, the cal- 
culation is repeated until the solution is converged up to 10"^ eV/ A'^. The lowest eigenvalue 
Ui is on average obtained after 1.5 iterations or equivalently 24 force evaluations. Compared 
to calculations based on the relative Lyapunov indicator, the present approach ensures sta- 
bility and high accuracy in the evaluation of ui which is essential in the estimation of the 
activation indicator. 
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The Monte Carlo method consists in samphng the trajectorial phase space endowed with 
parametrized probabihty density 



where fe, the reduced free energy of the biased trajectory ensemble, acts as a normalizing 
constant. The 6 parameter controls the biasing strength of the activation indicator in the 
path action 



Note that P^(2) = if xq ^ A : paths are all initiated from the reactant basin whatever 9. 
Besides, fe in f lT9|) depends not only on 9 but also on A and the inverse temperature /3. In 
the biased path ensemble, increasing the value of 9 decreases {C{z))g the ensemble average of 
the activation indicator. It thus enhances the occurrence of excursions into unstable regions 
of the potential energy landscape and the occurrence of active trajectories within the basin 
and the reactive trajectories crossing a saddle region. Conversely, decreasing the value 9 
favors the occurrence of inactive trajectories. 

The present ensemble is similar to the s-ensemble used by Chandler and co-workers'^^ to 
bias the occurence of dynamical heterogeneities in glasses. Their biasing field s couples to a 
dynamical activity defined to be an excitation indicator (number of spin flips in a trajectory) 
in the same way as the conjugate parameter 9 couples to C. 

The distribution is approximated by a Markov chain of M steps constructed using a 
Metropolis Monte Carlo method. As in TPS or Ref. |3^, the scheme consists in repeatedly 
attempting shooting and shifting moves. 

A. Shooting moves 

A shooting move consists in performing the following operations: 

(i) select a state = (qj", pj)'^ in the current trajectory by drawing an integer i uniformly 



(ii) generate trial momenta pe from current momenta pe and define trial state x^ = 




(19) 




(20) 



in [0,L]; 
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(iii) construct the trial trajectory z by propagating two trajectory segments from X£, one 
backward of duration ir and the other one forward of duration T — ^r; 

(iv) compute the acceptance probabihty of z 

Pace \z^ z\= min {1, /ia(xo) exp \ue{z) - ; (21) 

(v) draw a random number Randl G (0, 1] and accept the trial trajectory 'z if Randl < 
Pace \^ ^ z] otherwise the trial trajectory is rejected and the current trajectory is 
duplicated in the Markov chain. 

As detailed below, acceptance probability (|2T1) formally corresponds to the traditional 
Metropolis acceptance rate 



\z z\ = mm 



where Pgen \z z\ is the probability to generate z from z and vice versa for Pgon [2; ^ i] . 
Consequently, the probability fluxes Pace \z ^ z] Pgcn \z ^ z] P^ [z] from z to z and 
Pace [-2 ^ ^ Pgen ^ 2] P^ \z\ from z to z are balanced, and, shooting moves sample the 
equihbrium distribution P^. 

To show that the formal Metropolis rate ( 122|) and the rate ( 12T|) used in practice are 
equivalent, let us detail the probability of generating the trial trajectory z from z 

Pgen[? ^ Z\= Pcond(^|Xi')p {Vl ^ Vl) l—TT- (2^) 

The quantities {L + 1)~^ and p (p^ ^ p^) are the probabilities to draw ^ in step (i) and 
to generate the trial momenta p from p in step (ii) using a simple procedure proposed by 
Stoltz'^^ based on an Ornstein-Uhlenbeck process on the momenta 



P£ = ep£ + Vl - e25p^ (24) 

Here e is a mixing parameter and each 5pi/ variate is drawn from a Gaussian distribution 
of variance mj//3. Thus, the probability of generating p^ from p^ reads 



P (Pe ^ Pe) = exp 



2m. fl - e^] 
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where = Y[i=i a/ /3/[2mi(l — e^)??]. This form ensures that the probabihty fluxes between 
the current and perturbed momenta are balanced 

p (p£ ^ p^) ^ exp [-l3JC{pi)] 
pipi^Pi) exp[-/3/C(p^)]' 

As a corollary of the symplecticity property satisfied by (Vm)'^, the Jacobian and Gram 
determinants are equal to one whatever < i < L: the dynamics preserves the volume 
of any comoving infinitesimal hyperparallelepiped. As a result, the current and trial path 
probabilities conditioned on and respectively, [see equation (JTj)], become 

Pcond(2:|x<?) = Pcond(2|xo), (26a) 
Pcond(?|x^) = Pcond(?|xo)- (26b) 

Combining fl23|) . f l25|) and f l26|) yields the following ratio 



Pgen {Z ^Z) 


Pcond 1 


^z|xo) exp 


[-/3H(x,)] 


Pgen (Z ^ Z) 


Pcond 1 


^z|xo) exp 


[-/3H(5,)] 



(27) 



where we use the fact that the potential energy of the trial and current states at the shoot- 
ing index are equal to each other. Inserting the ratio ( 127|) in the Metropolis acceptance 
probability ( l22l) and resorting to the time symmetry property of the Verlet map yields the 
explicit form ( 12T|) of (iv). The small numerical drift on the Hamiltonian resulting from the 
integration of the position- Verlet scheme from Xq to x^ (^{(xf) ^ 'H(xo)) and from x^ to Xq 
('H(xo) — ^(x^)) has not been included in the Metropolis acceptance rate. This approxima- 
tion should not affect the temperature of the initial equilibrium state as each new momenta 
are generated from the exact Maxwell-Boltzmann distribution (using an Ornstein-Uhlenbeck 
process). The drift essentially affects the exactitude of Hamiltonian dynamics and must be 
kept small by using small enough a time-step. 

B. Shifting moves 

A shifting move is built upon a multiple proposal procedure^ and consists in performing 
the following operations 

(i) draw a random integer v in interval [0, L]; 
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(ii) propagate a first trajectory segment backward from xq for v steps and a second trajec- 
tory segment forward from x^^ for L — v steps; join them to z to form C = {x„}o<n<2L, 
the joint trajectory where x„ now denotes shifted state x„_,y; 

(iii) for each possible trial trajectory = {xj+n}o<n<L contained in the joint trajectory (, 
compute the associated selection probability 

Psei(^.IC) = hAiSc,) exp [UeiO - uoilj)] (28) 
where the joint action UelQ is defined by 

L 

exp [-UeiO] = /iyi(xj) exp [-ue{zj)] . 

j=0 

(iv) draw a random number Rand2 G (0, 1] and evaluate the lowest integer i such that 
E5=oPfei(?.lC)>Rand2. 

(v) add Zi to the Markov chain. 

To show that the shifting moves sample the distribution P^, we will prove that the 
transition probability fluxes between z = z,y and the selected trajectory Zi are balanced. 

The conditional probability to construct joint trajectory C starting from the current 
trajectory Zj reads 

PcondlCl^jO = Pbwd(C|j,5i)Pfwd(C|j + L,Sij+L)/iL + 1) 
= Pcond(CIXj)/ [Pcond(i'j|Xj)(L + 1)] 

where (L + is the probability to choose a particular u. The probability of each trial 
path Zj is P^^{zj) with the characteristic function Ha acting upon the initial state Xj [see 
Eq. (HI]. 

The marginal probability associated with a joint trajectory ( reads 

Pmarg(C) = Pcond (CI ) P A (^j ) • 

i=0 

which, for the position- Verlet map, simplifies into 

PLrg(C) = Pcond(CIXo) exp [fe - UoiC)] ■ (29) 
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wherein the indicator function Ha is hidden in the joint action Ug{(). The probabihty^S of 
selecting zj knowing the "joint" trajectory C is a posterior hkehhood given by Bayes relation 

pL(^.IC) = ^^^^P^^. (30) 

marg J 

The ratios of conditional probabilities that appear when developing (1301) simplifies into 1 
for the position- Verlet scheme (see fl26l) ). As a result, the posterior probability simplifies 
to ([28]). 

The multi-proposal shifting move consists of generating the joint trajectory ( from the 
current path z using Pcond(CI-2) and then selecting the next path of the Markov chain among 
the Zj proposals included in ( using posterior probability Pg^j. These shifting moves leave 
the probability distribution invariant because they satisfy the detailed balance condition 

PLl(2^.IC)Pcond(CI?.)PA(?.) = PLl^.lOPcondlCIS'.OPA (?.) 

where Psei(2i/|C)Pcond(C|2j) is the probability to transit from 'zj to z, and, conversely, 
Psei(^ilC)Pcond(C|2iy) is the probability to transit from z to Zj. 

A Monte Carlo cycle consists in performing a shooting move followed by a shifting 
one. The constructed sample consists of the trajectories obtained after shifting and is de- 
noted by {2;^, ■ ■ ■ , z*", ■ ■ ■ , z^}. The associated sample of joint trajectories is denoted by 



VI. ESTIMATING TRAJECTORY OBSERVABLES 

In the following, we employ a second control parameter a which formally may adopt the 
full range of possible values of 6. The intent of the two parameters are distinct however. 
The parameter a is meant to indicate the target distribution, i.e. the thermodynamic state 
that we are interested in measuring. As such the meaningful value of a will strictly be in 
the context of time-correlation functions. Measurements performed at a 7^ nevertheless 
enable one to assess the correctness of the results in term of fluctuations. Unlike the static 
a parameter, the 9 parameter refers to the various distributions P^ that our Monte Carlo 
algorithm has sampled. 
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A. Standard estimator 



The ensemble average of a trajectory observable O 



0{z)n{z)Vz 



(31) 



with respect to distribution given a sample of M trajectories {z'^}i<rn<M distributed 
according to can be approximated using the following online estimator 



M 



m=l 



0(2") exp 




Mexp 


fe - ne(z™) 





(32) 



The reduced free energy estimates fa and fg act as normalizing constants and may be 
obtained, up to an additive constant, by solving the equation /^*^ [1] = 1. 



B. Waste-recycling 

The statistical variance of the standard estimator can be reduced if one recycles infor- 
mation relative to unselected trial moves.— >^ Here, we will include information about all 
the unselected shifting moves contained in the joint trajectory C, (see Sec. IV BP . The de- 



sired waste- recycling estimator of trajectory observable O is constructed upon an ensemble 
average that is considered with respect to the marginal distribution 

{0)a = I E^(^^)Ps"i(^.IC)PL.g(C)^^C. (33) 

Resorting to Bayes relation (130|) . one checks that (133!) simplifies into the standard ensemble 
average fl3Tl) . Owing to this statistical equivalence, information from unselected trajectories 
can be retrieved^'29 provided that the sampler leaves a marginal probability distribution 
P^a^j.g invariant (with 6^ 7^ a in the general case). 

This requirement is met because our algorithm obeys an extended detailed balance be- 
tween any consecutive joint paths ( and The fulfilled detailed balance equation writes 

Pco„d(C'l4')vr [K> ^ Zu] Pfel(?.IC)PLrg(C) = Pcond(CF.)vr ^ 4'] Pfel (4 1 C')PL.g (O , 

(34) 
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where ly and i^, are the selected trajectories belonging to joint paths C, and C,' respectively (in 
the shifting). The quantity tt \z'^, -v- z^] denotes the transition probabilities (in the shooting) 
from Zy to i^/ and vice versa for vr [zy ^ zl,]. They obey detailed balance with respect to 
density. One checks that the detailed balance equation flMl) is indeed satisfied by (i) 
resorting to the shooting balance equation, (ii) inserting Bayes relation (!30l) within P^^^{zy\Q 
and Pfgi(i[,|C) and (iii) simplifying. Consequently, the distribution P^^rg l^ft invariant by 
our algorithm. 

Besides, the trajectory observable at Markov chain step m over the L + 1 proposals 
contained in joint trajectories C*" is averaged using the posterior probability P^^j of the 
target ensemble 



o: = 5^o(if)p^,i(ifin. 

j=0 

The ensemble average of a trajectory observable O with respect to marginal distribution 
-^marg ^xp (/„ — Ua) giveu a Sample of M joint trajectories distributed according to the 
marginal probability -P^arg ^^P if 9 ~ ^e) eventually be approximated using the fol- 
lowing waste-recycling estimator 



m=l 



exp 


'L 


-?7,(r) 


Mexp 


'fe 


-f/e(C'-) 



and again /„ — fg is approximated by solving J^g [1] 



C. Multi-state Bennett acceptance ratio 

In addition to the standard (|32|) and waste-recycling ( 135|) estimators, we also imple- 
ment the multi-state Bennett acceptance ratio (MBAR) method developed by Shirts and 
Chodera.— This method is based on Bennett acceptance ratio^ and extended bridge sam- 
pling^ techniques, and aims at minimizing the statistical variances associated with a series 
of simulations performed with K + 1 distinct values {0{k)}o<k<K of the control parameter 
6. The A;-th simulation provides one with a Markov chain consisting of Mg(fc) trajectories. 
Pooling all the data O"^ of the observable into a single chain of size M = ^e{k)^ the 
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waste-recycling estimate of observable O becomes 

M 



m=l 



exp 




-f/a(r) 




2^k= 


M0(fc) exp 





(36) 



A noticeable feature of MBAR is that the Markov chain origin of sample Z^fi IS cLIl irrelevant 
information. In practice, the free energy estimates { fe{k)}i<k<K are the solutions of the set 
of non-linear equations 

i?,Y,)(l) = l ^<k<K. 

The normalizing constants C6i(fc) = exp are evaluated up to a common multiplicative con- 
stant using the solver provided online.— In the waste- recycling context (Sec. IVIBp . the 
probability masses used in MBAR are the sampled marginal probabilities Pmarg(C'^) oc 
exp - Ue(^k){D] of m with Q<k<K. 

The state-to-state time-correlation functions C{nT) defined in Eq. [T]with t = nr and n < 
L can easily be computed from the MBAR estimator ( |36|) by setting a to 0. Noticing that 
Pseil^rlC"") simplifies into /^^(Xf)/ [Ei=i^A(X7)j and substituting /ib(5™+j) for ^(zf), 
the conditional time- correlation along joint trajectory ("^ reads 



(37) 



E-=o^A(5r) 

The posterior conditional average fl37|) characterizes the reactivity of the joint trajectory and 
is similar to the time average usually defining time-correlation functions.- Denoting B^^ (C) 
by C, the estimator of C{nr) reads 



C{nr) = J2 

m=l 



Without recychng the unselected candidates, the standard MBAR estimator would read 



{nr) exp 


7o- 






Ef=o exp 


f0{k) 


-Ueik)iC'^)_ 



M 



m=l 



0{z-^) exp 






Ef=o Me{k) exp 





(38) 



The standard MBAR form (138|1 differs from the waste-recycling form (136|) by the substitution 
of the probabihty ratios Pmarg/Pmarg = exp [fe — fa — Ue + Ua\ by the ratios P^/P^ = 
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exp [fe — fa — ue + Uq\, wherein uq = Ua + {0 — a)C Here, the fg(^k) estimates are solutions 
of Bf^^ll) = 1 and are determined up to a common additive constant. 

To our knowledge, this study is the first one to implement the MBAR reweighting scheme 
for waste-recycling. However, MBAR has already been applied to nonequilibrium path- 
ensemble averages^ for computing free energy differences or to Langevin dynamics for 
estimating time-correlation functions.— The path-sampling scheme of Sec. |V] can be im- 
plemented with Langevin dynamics^^ and even with stochastic dynamics based on discrete 
master equations.— The implementation with Langevin dynamics involves supplement- 
ing the Verlet splitting flTU]) with Ornstein-Uhlenbeck processes on the momenta^"— and 
generating the new proposals using Stoltz's Brownian tube technique.—'^ 



VII. APPLICATION TO VACANCY MIGRATION IN a-IRON 



SUNDAE is used to calculate the migration rate of a single vacancy in a-Fe, a crystalline 
phase of iron with body-centered cubic (BCC) structure. Atomic interactions of the model 
system are described by an embedded atom model potential.— The free energy profile asso- 
ciated with the jump of a neighboring atom into the vacancy has been calculated in Ref. |36 



using both Monte Carlo and the harmonic approximation for temperature ranging from 20 
K to 1000 K and with a lattice parameter (average first neighbor distance) a = 2.47A corre- 
sponding to the energy minimum at OK. The free energy profile exhibits a single free energy 
barrier for temperatures above T = 450 K, while for lower temperatures an intermediate 
metastable state appears, corresponding to the intra-site position of the jumping atom de- 



scribed in Ref. |53|. The overall process is clearly thermally activated for all temperatures. At 
500 K, the free energy barrier calculated in Monte Carlo is 0.501e\^. The reaction coordinate 
used to represent the migration was the projection of the atom jumping into the vacancy 
on an axis whose direction is parallel to the initial lattice sites of the vacancy and jumping 
atom (of the unrelaxed structure). 

The computational set-up in SUNDAE is as follows. Basin A and B are defined with 
respect to the underlying lattice whose sites are the atomic positions of the initial structure 
relaxed at K. The characteristic function Ha is 1 as long as all atoms are located within 
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a distance of 0.55A from their lattice site, otherwise it is 0. The characteristic function hs 
is 1 as soon as an atom is located beyond a distance of a/2 from its lattice site, otherwise 
it is 0. Trajectories are generated with time-step r = 2fs and consists of L = 300 steps. 
We have chosen the temperature of 500 K. We use K = 51 values of the 6 parameter and 
set d{k) = 6max X k/K with O^riax = 3.21. A preliminary simulation is done to equilibrate 
the trajectories at = ^^max- The final trajectory is active and is taken as starting point for 
the subsequent equilibration runs of 500 Monte Carlo cycles performed for all 6{k) values 
(0 < A; < K). A cycle is a shooting move followed by a shifting one. A few cycles are 
observed to inactivate the sampled trajectories at small 6 values. A production run consists 
of 2 X 10^ Monte Carlo cycles carried out for every 6 value independently: the sample sizes 
are thus Mg^k) = 2 x 10^ ioi < k < K . 

The value of the mixing parameters e (see Eq. [2^ is tuned to achieve the best trade-off 
between acceptance (which tends to 1 in the limit e = 1) and decorrelation (which arises in 
the limit e = 0). The increasing function e{6) = 1 — 2 x 100^^^^/^™=''= is a good compromise 
and ensures that the mean acceptance rate rj does not collapse below 10% with increasing 
6, as shown in Fig. [H 

Various MBAR estimates have been computed as a function of a. Free energy estimates 
fa are displayed in Fig. [2] together with their derivatives estimated from the ensemble average 
{C)a- The entropies Sa and their derivatives, displayed in Fig. [21 are estimated resorting 
to ensemble averages a{C)a — fa and —avaia i'^), respectively. The waste-recycling and 
standard forms agree with each other. 

The abrupt change in the slope of the reduced free energy, the inflection of the entropy 
and the sharp peak of the entropy derivative occurring at Oc = 1-48 are the signature of a 
phase transition in trajectory space.— The first phase consists of inactive and moderately 
active trajectories and is referred to as "the inactive phase". The second one, consisting 
of active trajectories exclusively, is favored as a increases and is referred to as the "active 
phase" . This behavior is emphasized by the histogram of the activation indicator displayed 
in Fig. m for three a values. The thermodynamic state is characterized by a single phase, 
inactive at a = and active at 2ac- In contrast at Oc, the inactive and active phases coexist. 
The phase coexistence around etc results in metastability issues in the simulations. 
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To show that the samples are indeed correlated, we have computed the effective number 
of uncorrected trajectories with respect to path observables £, Ce = d^^{Ue — Uq) and Ce- 
For any observable O and given the sample size Mg, this quantity is defined by^ 



E*^" ^ Mg — m Cov n(m) 



(39) 



^ Me Covo(O) 

m=l ^ ' 

where Covo is the autocovariance function of O along the stationary Markov chain: 

. Me 

Covo(m) = y rC"^C»"^+P - (^^o^+Pl . 

^ ' Mg-m ^ J 

p=m+l 

The sum in ( l39l) corresponds to the integrated fluctuation autocorrelation time. Figure [5] 
displays the effective number of uncorrelated trajectories and clearly shows the numerical 
slowing-down resulting from the phase transition. 

Let us now examine the standard errors for the online estimators associated with (£)q, 
which are displayed in Fig. |6l They were obtained by dividing each Markov chains into 200 
blocks consisting of 100 consecutive points, computing the statistical variance of the 200 
block estimates and eventually dividing their square roots by 200. Given the uncertainties, 
the online estimates are not reliable when a is close to ac- The same problematic trend 
holds for the estimates of the correlation observable l^^iz) = /iB(x(t)) displayed in Fig. [7] 
for t = T as a function of a. 

The standard error bars represented on the graphs in Figs. M and [7] for the two MBAR 
estimators have been obtained using the correlated data and thus underestimate the true un- 
certainty in MBAR. The use of uncorrelated subsamples in MBAR (as advocated in Ref. Iiol ) 
was problematic due to the small subsample sizes. Nevertheless, the MBAR uncertainties 
allow us to compare the efficiency of waste-recycling relatively to that of waste-disposal. We 
only observe a substantial improvement with waste-recycling at low a values. At a = 0, the 
standard error of B^'^iji^) is lower than the one of B^'^{h'^) by almost a factor 2, amounting 
to a simulation speed-up of almost 4. 

Since the average (/i^)o is the time-correlation function C(t), waste-recycling is useful 
for our problem. The information included in the unselected reactive trajectories is relevant 
and improves the accuracy of the estimates of C{t). This feature is well illustrated in 
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Fig. [HI the time-correlation function grows more linearly when estimated using the waste- 
recycling MBAR procedure than when estimated using the standard MBAR procedure. 
Concomitantly, the plateau value obtained for the derivative dC{t)/dt is considerably less 
noisy with waste-recycling. The plateau value corresponds to the phenomenological rate 
constant kA^B- 

We have also plotted in Fig. [8] the rates /ctst = A^'^ exp (— /JAF^i^b) given by 

transition state theory^*^ for comparison. The quantity A^'^ = 8 is the coordination number 
of BCC structure and Fa-^b is the free energy barrier from A to B, computed by Monte Carlo 
for a cell with 1023 atoms (MC) or calculated here using the classical harmonic 



in Ref. 
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approximation for the same cell of 127 atoms (HA). The small disagreements with the fcxsT 
value in the HA framework may originate from anharmonicities. The bigger disagreement 
with the A;tst value obtained by Monte Carlo may be explained by a size effect and by the 
fact that the employed reaction coordinate^ does not describe the reaction correctly, i.e. 
that the transmission factor^ associated with the reaction coordinate is small. Unfortunately, 



this transmission factor was not calculated in Ref. 
migration barrier exhibits a double hump. 



36 and is difficult to evaluate because the 



VIII. CONCLUSION 

Motivated by the ability of eigenvalue-foUowing^Si""— i^i^ and Lyapunov-weightingi^"— 
methods to locate saddle points in complex systems based on the topology of their energy 
surfaces, we developed and tested SUNDAE, a transition path sampling technique in which 
importance sampling in trajectory space does not confine the trajectory endpoints to a 
reactive state but rather uses the lowest eigenvalues of the Jacobian along the trajectory. 
In practice, the Lanczos algorithm is used to compute the lowest eigenvalue (as in the 
activation-relaxation technique^i2E) and the MBAR scheme is used in combination with 
waste-recycling to unbias the dynamically activated events and to estimate the state-to- 
state time- correlation function. 

The usefulness of the approach has been demonstrated on a practical example: the migra- 
tion of a vacancy in ct-Fe. The time-correlation function associated with a vacancy jump has 
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been computed. The plateau value of its time derivative has been found in good agreement 
with the rate constant given by the transition state theory. We also observe that the time 
derivative of the time-correlation function is considerably smoother when estimated with 
waste-recycling than without. Enabling MBAR to recycle the information contained in the 
unselected trajectories thus proved to be particularly relevant. 

To further speed-up the calculations, SUNDAE should be implemented in combination 
with a path replica exchange method^^ on the conjugate parameter controlling the 
mean value of the path action. This technique is routinely used to decorrelate samples in 
Markov chains^ and should alleviate the observed numerical slowing-down around the phase 
transition induced by the use of a conjugate parameter. The range of 9{k) values maximizing 
the efficiency of replica exchange simulations will likely have to span this transition. 

We believe that SUNDAE will find useful applications in a wide class of fields, spanning 
from molecular biophysics to physical metallurgy,— in which the knowledge of rate constants 
are important input parameters in kinetic Monte Carlo simulations or other meso-scale 
models. 

The authors are grateful to John D. Chodera (University of California, Berkeley) and 
Michael J. Shirts (University of Virginia) for their assistance with the implementations of 
MBAR and for stimulating comments. 
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FIG. 1. Mixing parameter e of Eq. [2H and mean acceptance rate as a function oi 9 = a. Each 77 
value is calculated from 2 x 10^ shooting attempts. At large 9 values, sampled paths are mainly 
active : if e is too small, then generated trial trajectories substancially decorrelate, become inactive 
and likely fail the Metropolis test. 
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FIG. 2. MBAR estimates of the reduced free energy and its derivative da fa as a function of 
a. The derivative is computed from the mean activation indicator (£)o. and denote the 
waste-recychng and standard MBAR estimators. 
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FIG. 3. MBAR estimates of the entropy Sa and of its derivative with respect to a. Same production 
run and same conditions as in Fig. [2j 
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FIG. 4. The quantity Probe (£) represents the probabihty to observe the value L during a simula- 
tion with 9 = a. The three histograms represents the distribution of the values of the activation 
indicator estimated using the standard MBAR estimator at a = 0, a = Oc, and a = 2qc, where 
Oc = 1.48. Same production run as in Fig. [2j 
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FIG. 5. Effective number of uncorrelated trajectories for all simulated values 0{k) of 6 with respect 
to the path functionals indicated in the legend. The ratio gg = Mq/Mq measures a statistical 
inefficiency {Mg = 2 x 10^ is the number of sampled trajectories for all the 0(A;)'s). 
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FIG. 6. Estimates of the mean activation indicator and of the standard error A{C)a as a function 
of a given by the MBAR estimators B^'^ and B^*^ and the onhne estimators /^^ and I^a- Same 
production run as in Fig. [2l 
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FIG. 7. Estimates of the time-correlation function {h^)a and of the standard error A{h^)a as a 
function of the unphysical a ensembles, using the MBAR estimators B'^'^ and B^'^ and the online 
estimators I^^^ and /^^a- Same production run as in Fig. [2j 
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FIG. 8. Standard and waste-recycling MBAR etimates of the correlation function C{t) = {hJ^)o 
(top panel) and of its time-derivative (bottom panel). TST rates are plotted for comparison with 
the obtained plateau value of the time derivative. 
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